#############################################
## PROGRAM NAME: 203a_res_msa.R            ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         201_af_msa_fin.Rda              ## 
##         202a_m1_msa.Rda -               ##
##         202a_m12_msa.Rda                ##
##                                         ##
## OUTPUTS:                                ##
##         Figures 4a-6a and 7             ##
##                                         ## 
## PURPOSE: Plot results MSA results       ##
##          (Figures 4a-6a and 7)          ## 
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

## load libraries ##
library(tidyverse) 
library(scales) 
library(ggtext) 

## compile the results ## 

ex.msm.res <- function(i, inlist) {
  
  res <- inlist[[i]][1]
  
  return(res)
  
}

ex.noadj.res <- function(i, inlist) {
  
  res <- inlist[[i]][3]
  
  return(res)
  
}

ex.adl.res <- function(i, inlist) {
  
  res <- inlist[[i]][4]
  
  return(res)
  
}

ex.fe.res <- function(i, inlist) {
  
  res <- inlist[[i]][5]
  
  return(res)
  
}

data_path <- # USER DEFINED PATH HERE #
setwd(data_path)

load("201_af_msa_fin.Rda")

## count non-missing obs ##

sum(!is.na(af.msa.fin$mgr_npov_2022))
sum(!is.na(af.msa.fin$mgr_hpov_2022))
sum(!is.na(af.msa.fin$mgr_hpov_cc_2022))

sum(!is.na(af.msa.fin$rb_npov_2022))
sum(!is.na(af.msa.fin$rb_hpov_2022))
sum(!is.na(af.msa.fin$rb_hpov_cc_2022))

sum(!is.na(af.msa.fin$mpv_npov_2022))
sum(!is.na(af.msa.fin$mpv_hpov_2022))
sum(!is.na(af.msa.fin$mpv_hpov_cc_2022))

sum(!is.na(af.msa.fin$eli_ushare_2022))
sum(!is.na(af.msa.fin$vli_ushare_2022))
sum(!is.na(af.msa.fin$li_ushare_2022))

## load results ## 

load("202a_m1_msa.Rda")
load("202a_m2_msa.Rda")
load("202a_m3_msa.Rda")
load("202a_m4_msa.Rda")
load("202a_m5_msa.Rda")
load("202a_m6_msa.Rda")
load("202a_m7_msa.Rda")
load("202a_m8_msa.Rda")
load("202a_m9_msa.Rda")
load("202a_m10_msa.Rda")
load("202a_m11_msa.Rda")
load("202a_m12_msa.Rda")


##
## msm results ## 
##

m1.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res1)

m1.msm.res.ov <- bind_rows(m1.msm.res) %>%
  filter(term == "thist")

m2.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res2)

m2.msm.res.ov <- bind_rows(m2.msm.res) %>%
  filter(term == "thist")

m3.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res3)

m3.msm.res.ov <- bind_rows(m3.msm.res) %>%
  filter(term == "thist")

m4.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res4)

m4.msm.res.ov <- bind_rows(m4.msm.res) %>%
  filter(term == "thist")

m5.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res5)

m5.msm.res.ov <- bind_rows(m5.msm.res) %>%
  filter(term == "thist")

m6.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res6)

m6.msm.res.ov <- bind_rows(m6.msm.res) %>%
  filter(term == "thist")

m7.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res7)

m7.msm.res.ov <- bind_rows(m7.msm.res) %>%
  filter(term == "thist")

m8.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res8)

m8.msm.res.ov <- bind_rows(m8.msm.res) %>%
  filter(term == "thist")

m9.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res9)

m9.msm.res.ov <- bind_rows(m9.msm.res) %>%
  filter(term == "thist")

m10.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res10)

m10.msm.res.ov <- bind_rows(m10.msm.res) %>%
  filter(term == "thist")

m11.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res11)

m11.msm.res.ov <- bind_rows(m11.msm.res) %>%
  filter(term == "thist")

m12.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res12)

m12.msm.res.ov <- bind_rows(m12.msm.res) %>%
  filter(term == "thist")

##
## noadj results ## 
##

m1.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res1)

m1.noadj.res.ov <- bind_rows(m1.noadj.res) %>%
  filter(term == "thist")

m2.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res2)

m2.noadj.res.ov <- bind_rows(m2.noadj.res) %>%
  filter(term == "thist")

m3.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res3)

m3.noadj.res.ov <- bind_rows(m3.noadj.res) %>%
  filter(term == "thist")

m4.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res4)

m4.noadj.res.ov <- bind_rows(m4.noadj.res) %>%
  filter(term == "thist")

m5.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res5)

m5.noadj.res.ov <- bind_rows(m5.noadj.res) %>%
  filter(term == "thist")

m6.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res6)

m6.noadj.res.ov <- bind_rows(m6.noadj.res) %>%
  filter(term == "thist")

m7.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res7)

m7.noadj.res.ov <- bind_rows(m7.noadj.res) %>%
  filter(term == "thist")

m8.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res8)

m8.noadj.res.ov <- bind_rows(m8.noadj.res) %>%
  filter(term == "thist")

m9.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = msa.res9)

m9.noadj.res.ov <- bind_rows(m9.noadj.res) %>%
  filter(term == "thist")

m10.noadj.res <- lapply(seq(1,100,by=1),
                        ex.noadj.res,
                        inlist = msa.res10)

m10.noadj.res.ov <- bind_rows(m10.noadj.res) %>%
  filter(term == "thist")

m11.noadj.res <- lapply(seq(1,100,by=1),
                        ex.noadj.res,
                        inlist = msa.res11)

m11.noadj.res.ov <- bind_rows(m11.noadj.res) %>%
  filter(term == "thist")

m12.noadj.res <- lapply(seq(1,100,by=1),
                        ex.noadj.res,
                        inlist = msa.res12)

m12.noadj.res.ov <- bind_rows(m12.noadj.res) %>%
  filter(term == "thist")


##
## adl results ## 
##

m1.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res1)

m1.adl.res.ov <- bind_rows(m1.adl.res) %>%
  filter(term == "zri_full_2022_st")

m2.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res2)

m2.adl.res.ov <- bind_rows(m2.adl.res) %>%
  filter(term == "zri_full_2022_st")

m3.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res3)

m3.adl.res.ov <- bind_rows(m3.adl.res) %>%
  filter(term == "zri_full_2022_st")

m4.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res4)

m4.adl.res.ov <- bind_rows(m4.adl.res) %>%
  filter(term == "zri_full_2022_st")

m5.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res5)

m5.adl.res.ov <- bind_rows(m5.adl.res) %>%
  filter(term == "zri_full_2022_st")

m6.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res6)

m6.adl.res.ov <- bind_rows(m6.adl.res) %>%
  filter(term == "zri_full_2022_st")

m7.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res7)

m7.adl.res.ov <- bind_rows(m7.adl.res) %>%
  filter(term == "zri_full_2022_st")

m8.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res8)

m8.adl.res.ov <- bind_rows(m8.adl.res) %>%
  filter(term == "zri_full_2022_st")

m9.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = msa.res9)

m9.adl.res.ov <- bind_rows(m9.adl.res) %>%
  filter(term == "zri_full_2022_st")

m10.adl.res <- lapply(seq(1,100,by=1),
                      ex.adl.res,
                      inlist = msa.res10)

m10.adl.res.ov <- bind_rows(m10.adl.res) %>%
  filter(term == "zri_full_2022_st")

m11.adl.res <- lapply(seq(1,100,by=1),
                      ex.adl.res,
                      inlist = msa.res11)

m11.adl.res.ov <- bind_rows(m11.adl.res) %>%
  filter(term == "zri_full_2022_st")

m12.adl.res <- lapply(seq(1,100,by=1),
                      ex.adl.res,
                      inlist = msa.res12)

m12.adl.res.ov <- bind_rows(m12.adl.res) %>%
  filter(term == "zri_full_2022_st")

## fixed effects ##

m1.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res1)

m1.fe.res.ov <- bind_rows(m1.fe.res) %>%
  filter(term == "zri_full_st")

m2.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res2)

m2.fe.res.ov <- bind_rows(m2.fe.res) %>%
  filter(term == "zri_full_st")

m3.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res3)

m3.fe.res.ov <- bind_rows(m3.fe.res) %>%
  filter(term == "zri_full_st")

m4.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res4)

m4.fe.res.ov <- bind_rows(m4.fe.res) %>%
  filter(term == "zri_full_st")

m5.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res5)

m5.fe.res.ov <- bind_rows(m5.fe.res) %>%
  filter(term == "zri_full_st")

m6.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res6)

m6.fe.res.ov <- bind_rows(m6.fe.res) %>%
  filter(term == "zri_full_st")

m7.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res7)

m7.fe.res.ov <- bind_rows(m7.fe.res) %>%
  filter(term == "zri_full_st")

m8.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res8)

m8.fe.res.ov <- bind_rows(m8.fe.res) %>%
  filter(term == "zri_full_st")

m9.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = msa.res9)

m9.fe.res.ov <- bind_rows(m9.fe.res) %>%
  filter(term == "zri_full_st")

m10.fe.res <- lapply(seq(1,100,by=1),
                     ex.fe.res,
                     inlist = msa.res10)

m10.fe.res.ov <- bind_rows(m10.fe.res) %>%
  filter(term == "zri_full_st")

m11.fe.res <- lapply(seq(1,100,by=1),
                     ex.fe.res,
                     inlist = msa.res11)

m11.fe.res.ov <- bind_rows(m11.fe.res) %>%
  filter(term == "zri_full_st")

m12.fe.res <- lapply(seq(1,100,by=1),
                     ex.fe.res,
                     inlist = msa.res12)

m12.fe.res.ov <- bind_rows(m12.fe.res) %>%
  filter(term == "zri_full_st")


## plot the results ##

msm.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)",
                                 "ELI supply",
                                 "VLI supply",
                                 "LI supply"),
                         coeff = c(mean(m1.msm.res.ov$estimate),
                                   mean(m2.msm.res.ov$estimate),
                                   mean(m3.msm.res.ov$estimate),
                                   mean(m4.msm.res.ov$estimate),
                                   mean(m5.msm.res.ov$estimate),
                                   mean(m6.msm.res.ov$estimate),
                                   mean(m7.msm.res.ov$estimate),
                                   mean(m8.msm.res.ov$estimate),
                                   mean(m9.msm.res.ov$estimate),
                                   mean(m10.msm.res.ov$estimate),
                                   mean(m11.msm.res.ov$estimate),
                                   mean(m12.msm.res.ov$estimate)),
                         lb = c(quantile(m1.msm.res.ov$estimate,0.025),
                                quantile(m2.msm.res.ov$estimate,0.025),
                                quantile(m3.msm.res.ov$estimate,0.025),
                                quantile(m4.msm.res.ov$estimate,0.025),
                                quantile(m5.msm.res.ov$estimate,0.025),
                                quantile(m6.msm.res.ov$estimate,0.025),
                                quantile(m7.msm.res.ov$estimate,0.025),
                                quantile(m8.msm.res.ov$estimate,0.025),
                                quantile(m9.msm.res.ov$estimate,0.025),
                                quantile(m10.msm.res.ov$estimate,0.025),
                                quantile(m11.msm.res.ov$estimate,0.025),
                                quantile(m12.msm.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.msm.res.ov$estimate,0.975),
                                quantile(m2.msm.res.ov$estimate,0.975),
                                quantile(m3.msm.res.ov$estimate,0.975),
                                quantile(m4.msm.res.ov$estimate,0.975),
                                quantile(m5.msm.res.ov$estimate,0.975),
                                quantile(m6.msm.res.ov$estimate,0.975),
                                quantile(m7.msm.res.ov$estimate,0.975),
                                quantile(m8.msm.res.ov$estimate,0.975),
                                quantile(m9.msm.res.ov$estimate,0.975),
                                quantile(m10.msm.res.ov$estimate,0.975),
                                quantile(m11.msm.res.ov$estimate,0.975),
                                quantile(m12.msm.res.ov$estimate,0.975)))


## no adj ##

noadj.res.in <- data.frame(var = c("Median gross rent (nat)",
                                   "Median gross rent (hpov)",
                                   "Median gross rent (hpov cc)",
                                   "Rent burden (nat)",
                                   "Rent burden (hpov)",
                                   "Rent burden (hpov cc)",
                                   "MPV (nat)",
                                   "MPV (hpov)",
                                   "MPV (hpov cc)",
                                   "ELI supply",
                                   "VLI supply",
                                   "LI supply"),
                           coeff = c(mean(m1.noadj.res.ov$estimate),
                                     mean(m2.noadj.res.ov$estimate),
                                     mean(m3.noadj.res.ov$estimate),
                                     mean(m4.noadj.res.ov$estimate),
                                     mean(m5.noadj.res.ov$estimate),
                                     mean(m6.noadj.res.ov$estimate),
                                     mean(m7.noadj.res.ov$estimate),
                                     mean(m8.noadj.res.ov$estimate),
                                     mean(m9.noadj.res.ov$estimate),
                                     mean(m10.noadj.res.ov$estimate),
                                     mean(m11.noadj.res.ov$estimate),
                                     mean(m12.noadj.res.ov$estimate)),
                           lb = c(quantile(m1.noadj.res.ov$estimate,0.025),
                                  quantile(m2.noadj.res.ov$estimate,0.025),
                                  quantile(m3.noadj.res.ov$estimate,0.025),
                                  quantile(m4.noadj.res.ov$estimate,0.025),
                                  quantile(m5.noadj.res.ov$estimate,0.025),
                                  quantile(m6.noadj.res.ov$estimate,0.025),
                                  quantile(m7.noadj.res.ov$estimate,0.025),
                                  quantile(m8.noadj.res.ov$estimate,0.025),
                                  quantile(m9.noadj.res.ov$estimate,0.025),
                                  quantile(m10.noadj.res.ov$estimate,0.025),
                                  quantile(m11.noadj.res.ov$estimate,0.025),
                                  quantile(m12.noadj.res.ov$estimate,0.025)),
                           ub = c(quantile(m1.noadj.res.ov$estimate,0.975),
                                  quantile(m2.noadj.res.ov$estimate,0.975),
                                  quantile(m3.noadj.res.ov$estimate,0.975),
                                  quantile(m4.noadj.res.ov$estimate,0.975),
                                  quantile(m5.noadj.res.ov$estimate,0.975),
                                  quantile(m6.noadj.res.ov$estimate,0.975),
                                  quantile(m7.noadj.res.ov$estimate,0.975),
                                  quantile(m8.noadj.res.ov$estimate,0.975),
                                  quantile(m9.noadj.res.ov$estimate,0.975),
                                  quantile(m10.noadj.res.ov$estimate,0.975),
                                  quantile(m11.noadj.res.ov$estimate,0.975),
                                  quantile(m12.noadj.res.ov$estimate,0.975)))

## adl ##

adl.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)",
                                 "ELI supply",
                                 "VLI supply",
                                 "LI supply"),
                         coeff = c(mean(m1.adl.res.ov$estimate),
                                   mean(m2.adl.res.ov$estimate),
                                   mean(m3.adl.res.ov$estimate),
                                   mean(m4.adl.res.ov$estimate),
                                   mean(m5.adl.res.ov$estimate),
                                   mean(m6.adl.res.ov$estimate),
                                   mean(m7.adl.res.ov$estimate),
                                   mean(m8.adl.res.ov$estimate),
                                   mean(m9.adl.res.ov$estimate),
                                   mean(m10.adl.res.ov$estimate),
                                   mean(m11.adl.res.ov$estimate),
                                   mean(m12.adl.res.ov$estimate)),
                         lb = c(quantile(m1.adl.res.ov$estimate,0.025),
                                quantile(m2.adl.res.ov$estimate,0.025),
                                quantile(m3.adl.res.ov$estimate,0.025),
                                quantile(m4.adl.res.ov$estimate,0.025),
                                quantile(m5.adl.res.ov$estimate,0.025),
                                quantile(m6.adl.res.ov$estimate,0.025),
                                quantile(m7.adl.res.ov$estimate,0.025),
                                quantile(m8.adl.res.ov$estimate,0.025),
                                quantile(m9.adl.res.ov$estimate,0.025),
                                quantile(m10.adl.res.ov$estimate,0.025),
                                quantile(m11.adl.res.ov$estimate,0.025),
                                quantile(m12.adl.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.adl.res.ov$estimate,0.975),
                                quantile(m2.adl.res.ov$estimate,0.975),
                                quantile(m3.adl.res.ov$estimate,0.975),
                                quantile(m4.adl.res.ov$estimate,0.975),
                                quantile(m5.adl.res.ov$estimate,0.975),
                                quantile(m6.adl.res.ov$estimate,0.975),
                                quantile(m7.adl.res.ov$estimate,0.975),
                                quantile(m8.adl.res.ov$estimate,0.975),
                                quantile(m9.adl.res.ov$estimate,0.975),
                                quantile(m10.adl.res.ov$estimate,0.975),
                                quantile(m11.adl.res.ov$estimate,0.975),
                                quantile(m12.adl.res.ov$estimate,0.975)))

## fe ##

fe.res.in <- data.frame(var = c("Median gross rent (nat)",
                                "Median gross rent (hpov)",
                                "Median gross rent (hpov cc)",
                                "Rent burden (nat)",
                                "Rent burden (hpov)",
                                "Rent burden (hpov cc)",
                                "MPV (nat)",
                                "MPV (hpov)",
                                "MPV (hpov cc)",
                                "ELI supply",
                                "VLI supply",
                                "LI supply"),
                        coeff = c(mean(m1.fe.res.ov$estimate),
                                  mean(m2.fe.res.ov$estimate),
                                  mean(m3.fe.res.ov$estimate),
                                  mean(m4.fe.res.ov$estimate),
                                  mean(m5.fe.res.ov$estimate),
                                  mean(m6.fe.res.ov$estimate),
                                  mean(m7.fe.res.ov$estimate),
                                  mean(m8.fe.res.ov$estimate),
                                  mean(m9.fe.res.ov$estimate),
                                  mean(m10.fe.res.ov$estimate),
                                  mean(m11.fe.res.ov$estimate),
                                  mean(m12.fe.res.ov$estimate)),
                        lb = c(quantile(m1.fe.res.ov$estimate,0.025),
                               quantile(m2.fe.res.ov$estimate,0.025),
                               quantile(m3.fe.res.ov$estimate,0.025),
                               quantile(m4.fe.res.ov$estimate,0.025),
                               quantile(m5.fe.res.ov$estimate,0.025),
                               quantile(m6.fe.res.ov$estimate,0.025),
                               quantile(m7.fe.res.ov$estimate,0.025),
                               quantile(m8.fe.res.ov$estimate,0.025),
                               quantile(m9.fe.res.ov$estimate,0.025),
                               quantile(m10.fe.res.ov$estimate,0.025),
                               quantile(m11.fe.res.ov$estimate,0.025),
                               quantile(m12.fe.res.ov$estimate,0.025)),
                        ub = c(quantile(m1.fe.res.ov$estimate,0.975),
                               quantile(m2.fe.res.ov$estimate,0.975),
                               quantile(m3.fe.res.ov$estimate,0.975),
                               quantile(m4.fe.res.ov$estimate,0.975),
                               quantile(m5.fe.res.ov$estimate,0.975),
                               quantile(m6.fe.res.ov$estimate,0.975),
                               quantile(m7.fe.res.ov$estimate,0.975),
                               quantile(m8.fe.res.ov$estimate,0.975),
                               quantile(m9.fe.res.ov$estimate,0.975),
                               quantile(m10.fe.res.ov$estimate,0.975),
                               quantile(m11.fe.res.ov$estimate,0.975),
                               quantile(m12.fe.res.ov$estimate,0.975)))


## set 1: median rents ## 

msm.fin.res1 <- msm.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 1 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 1 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 1 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res1 <- noadj.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 2 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 2 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 2 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res1 <- adl.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 3 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 3 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 3 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res1 <- fe.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 4 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 4 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 4 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res1 <- rbind(msm.fin.res1,
                  noadj.fin.res1,
                  adl.fin.res1,
                  fe.fin.res1)

fin.res1$Model <- factor(fin.res1$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

fin.res1$group <- factor(fin.res1$group,
                         levels = c("Above Average (N=349)",
                                    "High Poverty (N=221)",
                                    "High-Poverty Central City\n(N=207)")) 

levels(fin.res1$group) <- c(
  "<b>Above Average </b>(N=349)",
  "<b>High Poverty </b>(N=221)",
  "<b>High-Poverty<br>Central City </b>(N=207)"
)


## thanks to Google Gemini and Microsoft Copilot for en_dash formatting help below

en_dash_format <- function(x) {
  x <- scales::label_comma()(x)
  gsub("-", "\u2013", as.character(x))
}


ggplot(fin.res1, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = en_dash_format) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12),
    strip.text.y = element_text(size = 12, face = "bold")
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Monthly \nMedian Gross Rent	 ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))

## set 2: rent burden ##

msm.fin.res2 <- msm.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=349)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 1 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 1 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 1 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res2 <- noadj.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=349)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 2 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 2 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 2 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res2 <- adl.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=349)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 3 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 3 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 3 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res2 <- fe.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=349)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City\n(N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 4 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 4 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 4 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res2 <- rbind(msm.fin.res2,
                  noadj.fin.res2,
                  adl.fin.res2,
                  fe.fin.res2)

fin.res2$Model <- factor(fin.res2$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

fin.res2$group <- factor(fin.res2$group,
                         levels = c("Above Average (N=349)",
                                    "High Poverty (N=221)",
                                    "High-Poverty Central City\n(N=207)")) 

levels(fin.res2$group) <- c(
  "<b>Above Average </b>(N=349)",
  "<b>High Poverty </b>(N=221)",
  "<b>High-Poverty<br>Central City </b>(N=207)"
)


## thanks to @ Aaron for help removing leading zeros from y axis: https://stackoverflow.com/questions/42578262/remove-leading-zeros-in-ggplot2-scales

ggplot(fin.res2, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = c("0" = 0, ".01" = 0.01, ".02" = 0.02, ".03" = 0.03)) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    axis.ticks.x = element_blank(),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), # Adjust group label text
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Proportion \nRent-Burdened Households") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))



## set 3: MPV in impoverished neighborhoods ##

msm.fin.res3 <- msm.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=349)",
                           var == "MPV (hpov)" ~ "High Poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 1 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 1 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 1 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res3 <- noadj.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=349)",
                           var == "MPV (hpov)" ~ "High Poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 2 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 2 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 2 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res3 <- adl.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=349)",
                           var == "MPV (hpov)" ~ "High Poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 3 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 3 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 3 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.res3 <- fe.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=349)",
                           var == "MPV (hpov)" ~ "High Poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 4 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 4 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 4 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res3 <- rbind(msm.fin.res3,
                  noadj.fin.res3,
                  adl.fin.res3,
                  fe.fin.res3)

fin.res3$Model <- factor(fin.res3$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

fin.res3$group <- factor(fin.res3$group,
                         levels = c("Above Average (N=349)",
                                    "High Poverty (N=218)",
                                    "High-Poverty Central City (N=204)"))

levels(fin.res3$group) <- c(
  "<b>Above Average </b>(N=349)",
  "<b>High Poverty </b>(N=218)",
  "<b>High-Poverty<br>Central City </b>(N=204)"
)


ggplot(fin.res3, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = en_dash_format) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), # Adjust group label text
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Median \n Property Value ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## set 4: Supply ##

msm.fin.res4 <- msm.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 1 (N=349)",
                         var == "VLI supply" ~ "VLI supply 1 (N=349)",
                         var == "LI supply" ~ "LI supply 1 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res4 <- noadj.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 2 (N=349)",
                         var == "VLI supply" ~ "VLI supply 2 (N=349)",
                         var == "LI supply" ~ "LI supply 2 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res4 <- adl.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 3 (N=349)",
                         var == "VLI supply" ~ "VLI supply 3 (N=349)",
                         var == "LI supply" ~ "LI supply 3 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.res4 <- fe.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 4 (N=349)",
                         var == "VLI supply" ~ "VLI supply 4 (N=349)",
                         var == "LI supply" ~ "LI supply 4 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "FE")


fin.res4 <- rbind(msm.fin.res4,
                  noadj.fin.res4,
                  adl.fin.res4,
                  fe.fin.res4)

fin.res4$group <- factor(fin.res4$group,
                         levels = c("ELI supply (N=349)",
                                    "VLI supply (N=349)",
                                    "LI supply (N=349)"))

fin.res4$Model <- factor(fin.res4$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


levels(fin.res4$group) <- c(
  "<b>ELI supply </b>(N=349)",
  "<b>VLI supply </b>(N=349)",
  "<b>LI supply </b>(N=349)"
)


## thanks to Microsoft Copilot for the help with this formatting function 

en_dash_no_leading0 <- function(x, digits = 2) {
  # Base numeric formatting with consistent decimals
  lab <- scales::label_number(accuracy = 1 / (10^digits), trim = TRUE)(x)
  
  # Keep 0 as "0"
  lab[x == 0] <- "0"
  
  # Remove leading 0 for positive decimals: 0.xx -> .xx
  lab <- sub("^0\\.", ".", lab)
  # Remove leading 0 for negative decimals: -0.xx -> –.xx
  lab <- sub("^-0\\.", "\u2013.", lab)
  
  # For negative integers/other negatives: leading "-" -> en dash
  lab <- sub("^-", "\u2013", lab)
  
  lab
}



ggplot(fin.res4, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = c(0, -0.02, -0.04, -0.06, -0.08),
                     labels = en_dash_no_leading0) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), # Adjust group label text
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Proportion \n Housing Units") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))



## END OF PROGRAM ##

#sink()
